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1  Norms 


In  this  section  we  define  a  measure  of  mixing  that  does  not  necessarily  require  diffusion  to 
measure  the  amount  of  homogenization  that  occurs  during  the  mixing  process.  Recall  the 
advection-diffusion  equation 

c)0 

—  +  u-V9  =  kV29,  (1) 

where  6  is  a  concentration  field  in  a  finite  domain  U,  with  no-net-flux  boundary  conditions. 
We  assume  without  loss  of  generality  that 


9  dfi  =  0, 


and  define  the  L2-norm,  or  variance,  as 


(2) 


\\e\\l=  [  e2dn.  (3) 

Jn 

Recall  from  Lecture  1  that  the  variance  evolves  according  to 

±\\9\\l  =  -2K\\ve\\l  (4) 

and  decays  in  time  as  the  system  mixes.  The  variance  indicates  the  extent  to  which  the 
concentration  has  homogenized  and  is  thus  a  good  measure  of  the  amount  of  mixing  that 
has  occurred.  However,  the  variance  requires  knowledge  of  small  scales  in  6,  which  we  are 
not  necessarily  interested  in.  A  measure  of  how  well-mixed  the  concentration  is  does  not 
necessarily  require  knowledge  of  how  much  homogenization  has  occurred  due  to  diffusion  at 
small  scales.  This  is  more  in  keeping  with  the  definition  of  mixing  in  the  sense  of  ergodic 
theory  [2].  In  this  regard,  we  proceed  to  consider  the  pure  advection  equation 

^  +  u-  V6  =  0 
ot 

Note  that  in  this  case  equation  (4)  predicts  that  the  variance  satisfies 

>  11  =  0.  «» 


(5) 
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and  cannot  therefore  be  used  as  a  measure  of  mixing. 

The  advection  equation  (5)  takes  us  closer  to  the  ergodic  sense  of  mixing  in  which  we 
think  of  the  advection  due  to  the  velocity  field  as  a  time-dependent  operator  /S'*  :  T2  — >•  0 
that  moves  an  initial  patch  of  dye  according  to 

6q(x)  i->-  9(x,t)  =  St6o(x).  (7) 


If  we  consider  a  region  A  of  uniform  concentration  defined  by 


9q{x) 


l  if  x  e  A, 

0  otherwise, 


(8) 


then  the  volume  of  the  patch 

Vo\[6(x,t)]  =  Vo\(A),  (9) 

remains  constant  in  time  by  incompressibility.  We  can  associate  the  volume  of  the  patch 
with  the  Lebesgue  measure  and,  because  of  the  result  (9)  above,  St  is  measure-preserving. 
We  define  mixing  in  the  sense  of  ergodic  theory  by 

lim  Vol[A  n  S\B)}  =  Vol(A)Vol(R),  (10) 

f->-oo 

for  all  patches  A,  B  £  f l.  This  definition  follows  our  intuition  for  what  good  mixing  is. 
Referring  to  figure  1,  when  the  system  is  well-mixed  the  intersection  of  A  and  StB  is 
proportional  to  both  Vol(-A)  and  Vol(R).  Thus,  if  the  condition  (10)  holds  then  St  must 
spread  any  initial  patch  throughout  the  domain.  This  condition  is  referred  to  as  strong 
mixing  and  can  be  shown  to  imply  ergodicity. 

The  intersection  of  the  advected  patch  B  with  the  reference  patch  A  is  analogous  to 
projection  onto  L 2  functions.  This  motivates  the  following  weak  convergence  condition 


lim  (9(x,t),g)  =  0, 


(11) 


SlB 


Figure  1:  An  advected  patch  StB  that  has  undergone  strong  mixing.  At  late  times  the 
patch  covers  an  arbitrary  reference  patch  A. 
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for  all  functions  g  E  L2(D),  where  the  inner  product  is  defined  by 

(f,9)=  t  f(x)g(x)dn,  (12) 

Jn 

and  /  E  L2(Q )  if  fQ  |/|2  dfl  <  oo.  Weak  convergence  is  equivalent  to  mixing  as  a  conse¬ 
quence  of  the  Riemann-Lebesgue  lemma.  The  equivalent  conditions  (10)  and  (11)  require 
computing  over  all  patches  A  or  functions  g,  respectively.  Thus,  neither  of  these  conditions 
is  useful  in  practice.  However,  we  proceed  to  describe  a  theorem  that  shows  there  is  a 
simpler  way  to  determine  whether  or  not  weak  convergence  is  satisfied. 

Mathew,  Mezic  and  Petzold  [5]  introduced  the  mix-norm,  which  for  mean-zero  functions 
is  equivalent  to 

ll^llif-,/2  :=  1 1 V  1//26?  1 1 2  •  (13) 

Doering  and  Thiffeault  [1]  and  Lin,  Thiffeault  and  Doering  [3]  generalized  the  mix-norm  to 

P\\m--=  l|V90||2,  Q<  0,  (14) 

which  is  a  negative  homogeneous  Sobolev  norm.  This  norm  can  be  interpreted  for  negative 
q  via  eigenfunctions  of  the  Laplacian  operator.  For  example,  in  a  periodic  domain,  we  have 

ll«ll|,=£WW,  (15) 

k 

from  which  we  see  that,  for  q  <  0,  ||0||7.  smooths  6  before  taking  the  L2  norm.  The  theorem 
lim  ||0||jj9  =0,  q  <  0  <==>•  6  converges  weakly  to  0,  (16) 


Figure  2:  Comparison  of  the  mix- norms  for  a  flow  optimized  using  the  separate  methods  of 
optimal  control  and  optimal  instantaneous  decay.  Figure  from  Lin  et  al.  [3]. 
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t  =  0.6 


t  =  0.8 


Figure  3:  Evolution  of  the  concentration  field  for  the  flow  optimized  in  the  case  q  =  —  1  as 
computed  by  Lin  et  al.  [3]. 


due  to  Mathew,  Mezic  and  Petzold  [5]  and  Doering,  Lin  and  Thiffeault  [3]  shows  that  we 
can  track  any  mix- norm  to  determine  whether  a  system  is  mixing  (in  the  weak  sense). 
The  existence  of  this  quadratic  norm  makes  optimization  of  the  velocity  field  for  good 
mixing  considerably  easier.  Mathew,  Mezic,  Grivopoulos,  Vaidya  and  Petzold  [4]  have 
used  optimal  control  to  optimize  the  decay  of  the  q  =  —1/2  mix- norm.  Lin,  Doering  and 
Thiffeault  [3]  have  optimized  the  instantaneous  decay  rate  of  the  q  =  —  1  norm  using  the 
method  of  steepest  descent,  which  is  easier  to  compute  numerically  but  yields  suboptimal, 
but  nevertheless  very  effective,  stirring  velocity  fields.  A  comparison  of  the  methods  for 
optimized  mixing  is  shown  in  figure  2.  The  solid  line  decays  faster,  but  this  is  merely 
because  the  H cannot  be  compared  directly  with  FT-1/2.  The  corresponding  evolution  of 
the  concentration  field  for  the  case  q  =  —  1  from  Lin  et  al.  [3]  is  shown  in  figure  3. 
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2  Sources  and  Sinks 


Now  let  us  consider  the  situation  with  a  sources-sink  term  s(x,t), 

dt9  +  u-V9  =  nV2e  +  s{x,t),  (17) 

V  •  u  =  0. 

For  simplicity,  assume  that  f n  s(x,t)d£l  =  0.  Otherwise,  we  can  subtract  the  mean  of  9.  It 
is  convenient  to  think  of  sources  and  sinks  as  hot  and  cold  regions. 

Let  us  also  assume  that  our  sources  and  sinks  are  time- independent.  Then,  the  system 
eventually  achieves  a  steady  state  6(x)  that  satisfies 


u-V9  =  kV29  +  s. 

(18) 

We  define  the  operator 

C  =  y  •  V  —  kV2 

(19) 

so  that  (18)  can  be  written 

C9  =  s. 

(20) 

The  steady  solution  is  then 

9  =  C-ls 

(21) 

where  the  mean-zero  condition  on  9  makes  this  unique.  Note  that  k  /  0  is  needed  to 
achieve  a  steady  state.  So,  assuming  the  system  has  reached  a  steady-state,  we  have  to 
determine  how  we  measure  the  quality  of  mixing.  One  of  the  possible  ways  is  to  look  at 
the  norms  \\9 \\fjq,  where  q  =  0  represents  standard  derivation.  But  we  have  to  decide  what 
we  will  compare  to.  One  possibility  is  \\6 \\jjq/  |s \\jjq-  This  ratio  is  a  reasonable  choice,  but 
has  units  of  inverse  time.  It  is  preferable  to  use  a  dimensionless  quantity  for  measuring  the 
quality  of  mixing.  In  this  spirit,  we  define  mixing  enhancement  factors: 

\\d\\w 
q  P\W 

(22) 

where  9  is  the  purely-diffusive  solution  which  satisfies 

C9  =  s. 

Here,  £  =  — /-cV2  is  the  pure  diffusive  operator,  so  9  can  be  interpreted  as  the  solution  in  the 
absence  of  stirring.  Since  ||0||^9  is  usually  decreased  by  stirring,  eq  measures  the  enhance¬ 
ment  over  the  pure-diffusion  state.  Several  properties  are  given  in  Doering  and  Thiffeault  [1], 
Shaw,  Thiffeault  and  Doering  [7],  and  Thiffeault  and  Pavliotis  [8].  We  interpret  a  large  eq 
as  ‘good  stirring,’  since  it  in  that  case  the  norm  is  decreased  by  stirring. 

A  natural  question  is  whether  eq  can  ever  be  less  than  unity,  that  is,  if  stirring  can  ever 
be  worse  than  not  stirring.  Let’s  consider 


l|V0||2 
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Here, 

0  =  C~ls  =  (— kV2)  1  s  =  — k-1V-2s  ^  V0  =  -k^V-'s 
Also,  from  CO  =  s,  we  can  multiply  0  on  both  sides  and  take  spatial  average  and  then  get 

(0C0)  =  ( sO ) , 

where  (•)  =  fn  ■  dfl.  We  expand  the  left-hand  side: 

(0C0)  =  (Ou-VO)  -  n(0V20) 

=  (V-  (u0'2/2))  -  k(0V20) 

=  -k(0V20)  =  K(\V0\2). 

As  for  the  right-hand  side, it  can  be  written  as 

(Os)  =  (ov  •  v_1s>  =  -  (vo  •  v_1s>  =  k  (vo  ■  vo} , 

where  we  used 

0  =  C~1s  =  (-kV2)_1  s  =  -/e_1V-2s 
<=►  VO  =  -K-'V-'s. 

Recall  that  (|V0|2)  =  ll^ll2^.  Therefore, 

11^11^1  =  (vO-Vd)  <  j|V0||2||V0||2  =  11011^11011^. 

We  conclude  that 

II^IIri  <  Mfr  ^  S1  <1.  (23) 

This  is  somewhat  counter-intuitive  because  gradients  are  usually  increased  by  stirring.  How¬ 
ever,  the  gradients  in  a  steady-state  have  been  affected  by  diffusion. 

What  about  eq  for  values  of  q  other  than  1?  We  tried  and  failed  to  prove  eq  <  1,  simply 

because  it  is  not  true.  Following  a  challenge  by  Charlie  Doering  at  a  workshop  at  the  IMA 

in  2010,  Jeff  Weiss  came  up  with  something  like: 

u  =  (2  sin  x  cos  2 y  ,  —  cos  x  sin  2 y) ,  (24a) 

s  =  (cos  a:  —  ^)  siny  .  (24b) 

This  velocity  field  manages  to  concentrate  the  source  and  sink  distribution  more  than  dif¬ 
fusion  alone.  Streamlines  of  u  and  level  sets  of  s  are  shown  in  figure  4.  In  this  example,  we 
could  get  £o  —  0.978  and  e_i  ~  0.945,  which  are  slightly  less  than  1.  It  is  an  open  problem 
to  characterize  such  ‘unmixing’  flows. 
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(a) 


(b) 


x  x 

Figure  4:  The  pattern  of  velocity  field  and  source  for  the  ‘unmixing'  flow  and  source  distri¬ 
bution  (24). 

3  Optimization 

We  defined  the  mixing  enhancement  factors  based  on  Sobolev  norms.  Large  mixing  en¬ 
hancement  factor  indicates  good  mixing  for  a  given  source  and  sink  pattern.  One  of  the 
relevant  questions  in  this  step  is  what  kinds  of  flow  give  the  largest  eq  given  source  and  sink 
distribution  s(x). 

Here  is  a  simple  but  surprising  example.  The  source  and  sink  distribution  is  given 
by  s(x)  =  sinx  with  periodic  boundary  conditions  on  6.  The  optimal  solution  for  this 
source  and  sink  distribution  is  u  =  Ux,  which  is  constant  flow  from  the  hot  region  to  the 
cold  region  [6-8]  (figure  5).  This  example  demonstrates  that,  with  body  sources,  the  best 


Figure  5:  The  optimal  velocity  held  (solid  arrows)  for  the  source  distribution  s(x)  =  sinx. 
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Figure  6:  Optimal  stirring  velocity  field  (solid  lines)  for  the  source  s(x)  =  sin  x  sin  y  (colored 
background),  for  q  =  0  (left)  and  q  =  —  1  (right).  The  optimal  velocity  is  the  same  in 
both  cases  because  the  source  is  an  eigenfunction  of  the  Laplacian.  (Matlab  programs 
example  (1)  and  example  (2)  in  the  Appendix.) 


stirring  often  has  more  to  do  with  transport  than  with  creation  of  small  scales. 

More  generally,  we  have  to  solve  the  optimization  problem  numerically.  Figure  6  shows 
contours  of  the  streamfunction  for  the  optimal  stirring  velocity  (lines)  for  a  source  s(x)  = 
sin  x  sin  y,  for  q  =  0  and  q  =  —  1.  The  optimal  velocity  fields  are  identical  for  the  two  values 
of  q,  because  the  source  is  an  eigenfunction  of  the  Laplacian. 

Contrast  this  to  the  optimal  solutions  in  figure  7,  for  the  source  distribution  s(»)  = 
cos  x  cos  y  +  cos  3 y  +  (1/4)  sin3y.  This  source  is  not  an  eigenfunction  of  the  Laplacian,  and 
we  expect  optimal  solutions  to  depend  on  q.  Comparing  the  left  (q  =  0)  and  right  ( q  =  —  1) 
figures,  we  see  this  is  indeed  the  case,  though  the  difference  in  this  case  is  fairly  small. 

Finally,  given  an  optimization  code,  it  is  simple  to  turn  it  around  to  anti-optimize,  that 
is,  find  the  worst  stirring  velocity  for  a  given  source  distribution.  Figure  8  shows  this  for 
the  source  (24b)  and  q  =  0.  Note  how  the  velocity  field  seems  to  work  to  concentrate  the 
source  sink,  thereby  increasing  the  variance.  The  efficiency  for  this  anti-optimal  solution 
is  £q  =  0.9736,  which  is  not  much  lower  than  Jeff  Weiss’s  unoptimized  flow  (24a),  which 
had  £o  —  0.978. 

To  reproduce  the  5  figures  in  this  section  run  the  program  example(n)  in  the  Appendix, 
where  n  is  a  number  from  1  to  5. 


Appendix:  Matlab  code 

1  Program  file  example. m 

function  example (ex) 

7.  util  folder  contains  Diffmat2.m,  fourdif.m,  refine.m,  refine2.m,  refinek.m 
7.  Contact  author  to  obtain  these  functions, 
addpath  util 

if  nargin  <  1 ,  ex  =  1 ;  end 
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Figure  7:  Optimal  stirring  velocity  field  (solid  lines)  for  the  source  cos  x  cos  y  +  cos  3 y  + 
(1/4)  sin3y  (colored  background),  for  q  =  0  (left)  and  q  =  —  1  (right).  The  optimal  velocities 
are  different  since  the  source  is  not  an  eigenfunction  of  the  Laplacian.  (Matlab  programs 
example  (3)  and  example  (4)  in  the  Appendix.) 


Figure  8:  Optimal  ‘unmixing’  solution  for  the  source  (cosx  —  sin  y,  with  mixing  effi¬ 
ciency  £q  =  0.9736.  (Matlab  program  example  (5)  in  the  Appendix.) 
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N  =  11;  °/0  Number  of  gridpoints 
L  =  1;  kl  =  2*pi/L;  #/«  Physical  size  of  domain 
x  =  L*(0:N-1)/N;  y  =  x’ ;  [xx,yy]  =  meshgrid(x,y) ; 

switch  ex 

•/.  Examples  1  and  2  use  the  same  cellular  source,  for  q=0  and  q=-l. 

°/0  Since  the  source  is  an  eigenfuntion  of  the  Laplacian,  the  optimal 
°/e  flow  is  the  same  for  q=-l. 
case  {1,2} 

src  =  cos (kl*xx) . *cos (kl*yy)  *  2/L; 

psiO  =  sin(kl*xx) . *sin(kl*yy)  *  l/sqrt(2)/pi; 

kappa  =  .1;  q  =  1-ex; 

l  Examples  3  and  4  use  the  same  two-mode  source,  for  q=0  and  q=-l. 

I  Since  the  source  is  not  an  eigenfuntion  of  the  Laplacian,  so  the  optimal 
1  flow  is  different  for  q=-l. 
case  {3,4} 

src  =  (cos(kl*xx) . *cos(kl*yy)+cos(3*kl*yy)+. 25*sin(3*kl*yy) )  *  4*sqrt (2) /5/L 
psiO  =  sin(kl*xx) . *sin(kl*yy)  *  l/sqrt(2)/pi; 
kappa  =  .1;  q  =  3-ex; 
case  5 

%  Unmixing  solution 

src  =  (cos(kl*xx)  -  . 5) . *sin(kl*yy)  *  2*sqrt (2/3)/L; 
psiO  =  sin(kl*xx) . *sin(2*kl*yy)  *  1/sqrt (5) /pi ; 

scalefac  =  -N~2;  °/0  Set  scale  factor  negative  to  minimize  instead 
kappa  =1/4;  q  =  0; 
end 

if  ~exist( ’scalefac’ ) ,  scalefac  =  N~2;  end 
[psi,Effq]  =  velopt (psiO, src, kappa, q,L, scalefac) ; 
fprintf  (1,  ’Ef  f  _°/od=0/0f  \n’  ,q,Effq) 
figure (1) 

Nplot  =  64;  l  Interpolate  solution  for  plotting 
psir  =  ref ine2(psi , Nplot) ;  srcr  =  ref ine2(src, Nplot) ; 
xplot  =  L* (0:Nplot-l) /Nplot ;  yplot  =  xplot’; 
imagesc (xplot ,yplot , srcr) ,  colorbar,  hold  on 
contour (xplot , yplot ,psir, 10, ’EdgeColor’ , ’k’ ) ,  hold  off 


2  Program  file  velopt. m 

function  [psi,Effq]  =  velopt (psiO , src , kappa, q,L, scalefac) 

°/«  Problem  parameters  for  Matlab’s  optimizer  fmincon. 
psiO  =  psi0(:);  problem. xO  =  psiO (2 rend); 

problem. objective  =  @(x)  normHq2(x, src, kappa, q,L, scalefac) ; 
problem. nonicon  =  @(x)  nonlcon(x, src, kappa, q,L, scalefac) ; 
problem. solver  =  ’fmincon’; 

problem. options  =  optimset ( ’Display ’,’ iter ’, ’TolFun’ , le-10 ,.. . 

’GradObj ’ , ’on’ , ’GradConstr’ , ’on’ , . . . 

’ algorithm ’ , ’ interior-point ’ ) ; 

[psi,Hq2]  =  fmincon (problem) ; 

1  Mixing  efficiency:  call  normHq2  with  no  flow  to  get  pure-conduction  solution 
Effq  =  sqrt (normHq2 (zeros (size (psi) ) , src, kappa, q,L, scalefac)  /  Hq2) ; 

psi  =  reshape  (  [0;  psi]  ,  size  (src))  ;  °/«  Convert  psi  back  into  a  square  grid 


%========================================================= 

function  [varargout]  =  normHq2 (psi , src , kappa, q,L , scalefac) 

N  =  size(src,l);  src  =  src(r); 
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7,  2D  Differentiation  matrices  and  negative-Laplacian 
[Dx,Dy ,Dxx,Dyy]  =  Dif fmat2(N,L) ;  mlap  =  -(Dxx+Dyy); 

if  q  ~=  0  &&  q  ~=  -1,  error(’This  code  only  supports  q  =  0  or  -I.’);  end 


psi  =  [0;psi] ;  ux  =  Dy*psi;  uy  =  -Dx*psi; 

ugradop  =  diag ( sparse (ux) ) *Dx  +  diag ( sparse (uy) ) *Dy; 


if  q  ==  0 

Aop2  =  (-ugradop  +  kappa*mlap) ; 
elseif  q  ==  -1 

Aop2  =  mlap* (-ugradop  +  kappa*mlap) ; 
end 

Aopl  =  (ugradop  +  kappa*mlap) *Aop2 ; 

7,  Solve  for  chi,  dropping  corner  point  to  fix  normalisation, 
chi  =  [0;  Aopl (2 : end, 2 : end)  \  src(2:end)]; 
theta  =  Aop2*chi; 

7,  The  squared  H~q  norm  of  theta. 

varargout{l}  =  L'‘2*sum (theta.  "2) /N~2  *  scalefac; 

if  nargout  >  1 

7.  Gradient  of  squared-norm  Hq2. 

gradHq2  =  2* ( (Dx*theta) . * (Dy*chi)  -  (Dy*theta) . * (Dx*chi) ) ; 
varargout{2}  =  gradHq2(2 : end)  /  N~2  *  scalefac; 
end 


%============================================================= 

function  [c , ceq,gc ,gceq]  =  nonlcon(psi , src , kappa, q,L, scalefac) 

psi  =  [0;psi] ;  N  =  size(src,l); 
c  =  []  ;  gc  =  []  ; 

[Dx,Dy ,Dxx,Dyy]  =  Dif fmat2(N,L)  ;  7«  2D  Differentiation  matrices 
U2  =  L~2* (sum( (Dx*psi) . ~2  +  (Dy*psi)  .  ~2) /N~2)  ; 
ceq(l)  =  (U2-1)  *  scalefac; 

if  nargout  >  2 

7.  Gradient  of  constraints 
mlappsi  =  - (Dxx+Dyy) *psi; 

gceq(:,l)  =  2*mlappsi (2 : end)  /  N~2  *  scalefac; 
end 
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